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SSUSI  Auroral  E-region  Algorithm 


1.  Introduction 

The  SSUSI  Auroral  E-region  algorithm  (SSUSIAE)  characterizes  electron  and  proton 
precipitation  in  the  diffuse  aurora  using  select  ultraviolet  (UV)  emissions  produced  in  the  E-region 
of  the  auroral  ionosphere.  SSUSIAE  also  calculates  N^E  (the  peak  electron  density  of  the  E- 
region)  and  h^E  (the  height  of  the  E-region)  based  on  the  deduced  particle  precipitation 
characteristics. 

We  have  continued  to  base  the  algorithm  exclusively  on  Maxwellian  incident  spectra 
(appropriate  to  the  continuous  aurora)  rather  than  using  different  spectra  for  different  kinds  of 
aurora.  We  are  aware  that  discrete  auroral  arcs  are  produced  by  particles  with  a  narrower 
spectrum  (usually  represented  by  Gaussians),  but  we  feel  that  more  work  needs  to  be  done  on  the 
problem  of  automatic  discrimination  of  discrete  arcs  from  the  continuous  aurora.  We  anticipate 
that  some  sort  of  discrimination  will  be  included  in  the  multi-sensor  version  of  this  algorithm. 

We  have  continued  to  ignore  the  01  1356  A  emission  in  the  auroral  algorithm.  We  feel 
that  its  usefulness  as  an  operational  tool  in  the  single  sensor  algorithm  is  quite  limited.  However, 
we  expect  to  make  use  of  it  in  the  multi-sensor  version  of  the  algorithm. 


2.  Algorithm  Description 


SSUSIAE  is  designed  to  work  with  UV  images.  However,  except  for  the  geocoronal 
background  subtraction  required  for  the  Lyman  a  image,  and  the  dayglow  background 
subtraction  required  for  the  dayside  portion  of  the  auroral  oval,  it  handles  each  image  “pixel  by 
pixel.”  From  the  three  UV  images  that  it  accepts  as  input,  the  algorithm  produces  four  images  of 
particle  precipitation  parameters  and  two  “images”  of  vertical  electron  density  profile  parameters. 
As  delivered,  the  algorithm  assumes  that  the  three  images  are  contained  in  a  single  file  in  which 
each  pixel  is  tagged  by  the  UT  of  the  measurement,  the  geographic  and  geomagnetic  location,  and 
the  intensities  of  the  three  emission  features  required  by  the  algorithm.  The  input  routine  is  a 
separate  subprogram  and  can  be  easily  modified  to  accommodate  any  alternative  file  format. 

The  three  UV  emission  features  used  by  SSUSIAE  are: 

1.  =  Lyman  a  (1216  A) 

2.  LBHl  =  LBH  1450  A  band  (1400-1500  A) 

3.  LBH2  =  LBH  1725  A  band  (1650-1800  A) 

For  each  pixel,  SSUSIAE  determines  parameters  that  describe  the  incident  spectra  of  precipitating 
electrons  and  protons.  Since  the  purpose  of  this  algorithm  is  to  predict  behavior  in  the  diffuse 
auroral  E-region,  the  incident  particle  energy  spectra  are  assumed  to  have  a  Maxwellian  form  that 
is  described  by  two  parameters:  the  characteristic  energy  and  the  energy  flux  Q.  Therefore, 
SSUSIAE  calculates  four  precipitation  parameters  per  pixel: 

1.  £o,e.  the  characteristic  energy  of  the  incident  electrons 

2.  Q^,  the  energy  flux  of  the  incident  electrons 

3.  £■(,  p,  the  characteristic  energy  of  the  incident  protons 

4.  0p,  the  energy  flux  of  the  incident  protons 

Additionally,  SSUSIAE  determines  the  ionospheric  parameters  N^E  and  h^E  for  each 
pixel. 

The  algorithm  also  estimates  the  uncertainty  in  each  of  the  output  parameters  due 
uncertainties  of  the  input  data  and  the  internal  uncertainties  of  the  algorithm  itself  Many 
analytic  functions  that  form  the  auroral  algorithm  are  in  the  form  of  power  series: 
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where  /  is  either  the  yield  of  an  emission  feature  or  the  characteristic  energy  of  a  precipitating 
particle  spectrum,  a  is  either  an  intensity  ratio  or  a  characteristic  energy.  Generally  the 
characteristic  energy  is  first  determined  from  an  intensity  ratio.  Then  the  yield  of  an  emission 
feature  is  determined  from  the  characteristic  energy.  The  yield  is  used  to  determine  the  energy 
flux.  The  yield  and  energy  flux  may  be  used  to  subtract  the  contribution  of  one  kind  of  particle 
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(e.g.  protons)  from  the  total  intensity  to  reveal  the  contribution  by  another  kind  of  particle  (e.g., 
electrons).  So  the  algorithm  relies  on  repeated  evaluations  of  several  formulae  like  Equation  (1). 
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where  o  is  the  covariance  of  c„  andc^ .  L  =  ®  L  the  covariance  of  c„  and  a,  and  o  is  the 
variance  of  a.  By  definition,  a  vanishes,  so  Equation  (2)  collapses  to 
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The  are  provided  by  the  least  squares  fitting  process  that  produces  the  coefficients  c„  of 
Equation  (1).  When  a  is  an  intensity  or  intensity  ratio  or  other  measured  quantity,  a^^is  simply 
the  measurement  uncertainty.  When  a  is  a  derived  quantity,  its  uncertainty  is  itself  estimated 
using  an  equation  of  the  form  of  Equation  (3).  Finally,  the  uncertainty  in /itself  (rather  than  In  f) 
is  calculated  from 


«^/  =  ^in(/)/  (4) 

Note  that  the  accuracy  and  validity  of  the  above  uncertainty  estimates  depend  on  the 
manner  in  which  the  series  are  determined.  The  actual  intensities  depend  on  a  number  of 
parameters  that  are  not  known,  primarily  the  precise  composition  of  the  neutral  atmosphere. 
Although  model  atmospheres  can  give  reasonable  estimates  of  the  neutral  composition,  the  aurora 
itself  modifies  the  composition,  which  limits  the  level  of  accuracy  achievable  by  a  model 
atmosphere.  If  the  covariances  a  are  to  be  useful  as  an  estimate  of  the  uncertainty  in  the 
coefficients  due  to  the  variability  of  “uncontrolled”  parameters,  the  simulated  data  used  to 
produce  Equation  (1)  must  span  the  natural  range  of  these  parameters  but  not  give  too  much 
weight  to  extreme  values.  Realistically,  the  uncertainty  estimates  must  be  considered  crude 
estimates  and  should  not  be  over-interpreted. 

The  top  level  flow  diagram  for  SSUSIAE  appears  in  Figure  1.  The  flow  diagram  for  the 
main  part  of  the  algorithm,  AURORA,  appears  in  Figure  2.  The  sections  that  follow  describe 
each  of  the  major  components  of  SSUSIAE  in  some  detail. 

It  should  be  noted  that  this  algorithm  was  implemented  in  FORTRAN  with  the 
understanding  that  it  would  be  translated  into  ADA.  Therefore,  the  coding  was  done  to  make  the 
data  flow  explicit  and  easily  understandable,  rather  than  to  take  fullest  advantage  of  specific 
FORTRAN  structures. 
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It  should  be  further  noted  that  this  algorithm  assumes  that  the  image  has  been  smoothed, 
geolocated,  and  rectified.  That  is,  it  assumes  that 

1.  any  averaging  or  summing  required  to  reduce  the  statistical  noise  to  an  acceptable 
level  has  been  done, 

2.  each  pixel  is  tagged  with  the  correct  geographic  and  geomagnetic  coordinates, 
including  solar  zenith  angle, 

3.  region  boundaries  have  been  determined  so  that  each  pixel  is  also  tagged  with  the 
region  in  which  it  lies,  and 

4.  the  intensities  have  been  converted  from  slant  intensities  to  equivalent  vertical 
intensities. 
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SSUSIAE 

SSUSI  Auroral  E-layer  Algorithm 
Top  Level  Flow  Diagram 


READGP: 
read  geophysical 
parameters 


disk  Tile: 
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parameters 


READIMG: 

read  SSUSI  image  files 


disk  file: 
SSUSI  images 


_ _ t 

GEOCRNA: 
determine  geocoronal 
background 


DAYGLO: 

ij  determine  the  dayglow  | 
background 


_ ♦ 

AURORA: 

executive  subprogram 
of  auroral  algorithm 
(see  Fig.  2) 


Figure  1.  The  top  level  flow  diagram  for  SSUSIAE. 
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Figure  2.  Flow  diagram  for  AURORA  and  its  subsidiary  subroutines. 


6 


2.1  COMMON /IMAGE/ 


Because  the  image  data  is  used  unaltered  by  a  number  of  subroutines,  it  has  been  placed  in 
the  block  COMMON  /IMAGEy.  The  contents  of  /IMAGE/  conform  as  closely  as  possible  to  the 
expected  contents  of  the  SSUSI  Image  Mode  EDR  Data  Elements.  However,  /IMAGE/ 
represents  a  subset  of  the  SSUSI  EDR.  The  following  variables  and  arrays  are  contained  in 
/IMAGE/: 


DAY; 

UT; 

NPIXEL: 

NROW; 

NCOL: 

ROW; 

COLUMN: 

GLAT: 

GLON: 

MLAT: 

MLON; 

SZA; 

SLT: 

REGION: 

11216: 

11450: 

11725; 

C1216; 

C1450: 

C1725: 

U1216: 

U1450: 

U1725: 


The  day  of  the  year  of  the  image  [INTEGER,  unitless] 

The  universal  time  of  the  image  (presumably  the  UT  of  the  center  of  the  image) 
[REAL,  hours] 

The  actual  number  of  pixels  in  the  image.  [INTEGER,  unitless] 

The  number  of  rows  in  the  image.  [INTEGER,  unitless] 

The  number  of  columns  in  the  image.  [INTEGER,  unitless] 

The  row  number  of  each  pixel.  [INTEGER  ARRAY,  unitless] 

The  column  number  of  each  pixel.  [INTEGER  ARRAY,  unitless] 

The  geographic  (or  geodetic)  latitude  of  each  pixel.  [REAL  ARRAY,  degrees] 
The  geographic  (or  geodetic)  longitude  of  each  pixel.  [REAL  ARRAY, 
degrees] 

The  geomagnetic  latitude  of  each  pixel.  [REAL  ARRAY,  degrees] 

The  geomagnetic  longitude  of  each  pixel.  [REAL  ARRAY,  degrees] 

The  solar  zenith  angle  of  each  pixel.  [REAL  ARRAY,  degrees] 

The  solar  local  time  of  each  pixel.  [REAL  ARRAY,  hours] 

The  region  flag  for  each  pixel.  [INTEGER  ARRAY,  unitless] 

The  Lyman  a  intensity  for  each  pixel.  [REAL  ARRAY,  Rayleighs] 

The  LBH  1  (1400  -  1 500  A)  intensity  for  each  pixel.  [REAL  ARRAY, 
Rayleighs] 

The  LBH  2  (1650  -  1800  A)  intensity  for  each  pixel.  [REAL  ARRAY, 
Rayleighs] 

The  actual  number  of  counts  due  to  Lyman  a  for  each  pixel.  [REAL  ARRAY, 
counts] 

The  actual  number  of  counts  due  to  LBH  1  for  each  pixel.  [REAL  ARRAY, 
counts] 

The  actual  number  of  counts  due  to  LBH  2  for  each  pixel.  [REAL  ARRAY, 
counts] 

The  uncertainty  in  the  Lyman  a  intensity  for  each  pixel.  [REAL  ARRAY,  %] 
The  uncertainty  in  the  LBH  1  intensity  for  each  pixel.  [REAL  ARRAY,  %] 

The  uncertainty  in  the  LBH  2  intensity  for  each  pixel.  [REAL  ARRAY,  %] 


The  contents  of  /IMAGE/  are  read  from  a  disk  file  by  READIMG  (Section  2.5). 
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2.2  COMMON /ECOEFF/ and  COMMON /PCOEFF/ 


Because  the  coefficients  in  the  analytic  representations  of  yield  curves,  etc.,  that  make  up 
the  basic  algorithm  are  read  by  READCO  (Section  2.3)  and  used  only  in  the  subprogram 
PSEVAL,  which  is  deeply  nested,  they  are  also  stored  in  block  COMMON.  The  electron 
coefficients  are  stored  in  /ECOEFF/,  while  the  proton  coefficients  are  stored  in  /PCOEFF/.  The 
analytic  representations  and  coefficients  are  described  in  Section  2.3. 


2.3  SUBROUTINE  READCO 

SUBROUTINE  READCO  is  called  from  the  main  program  with  one  argument; 

LUN:  The  logical  unit  number  of  the  file  containing  the  data  that  is  to  be  read  (input) 

READCO  reads  the  coefficients  of  the  analytical  models  of  yield  vs.  characteristic  energy  and 
characteristic  energy  vs.  LBH  ratio.  These  coefficients  are  used  by  the  FUNCTION  subprogram 
PSEVAL  to  evaluate  the  yield  curves  and  other  functional  relationships  used  in  the  algorithm. 
The  coefficients  are  stored  in  a  common  block  rather  than  returned  as  an  argument. 

In  the  current  version  of  the  algorithm,  the  analytical  models  all  have  the  same  form: 

ln[/(a)]  =  (5) 

rt=0 

When / represents  Zi  yield  (the  optical  intensity  produced  by  an  energy  flux  of  lerg  cm'^  s  '),  the 
argument  a  is  the  characteristic  energy  of  the  precipitating  particles  producing  the  emission  (or  its 
inverse).  When / represents  the  characteristic  energy  of  the  precipitating  particles,  the  argument  a 
is  the  ratio  of  the  LBH  1400-1500  A  to  LBH  1650-1800  A  band  intensities  (or  its  inverse).  The 
identity  of  the  precipitating  particles  is  indicated  by  the  first  superscript  on  the  coefficients  c:  b  = 
e  for  electrons,  and  b  =  p  for  protons.  The  second  superscript  on  the  coefficients  is  merely  an 
integer  used  to  distinguish  the  coefficients  of  each  series.  For  protons: 

/  =  1  for  1^450,  the  yield  of  LBH  1400-1500  A 
/•  =  2  for};^7^^,  the  yield  of  LBH  1650-1800  A 
/  =  3  for}^2?g,  the  yield  of  Lyman  a  (1216  A) 

;■  =  4  for£o  p,  the  characteristic  energy  of  the  incident  proton  spectrum 
For  electrons: 

/•  =  1  forYl;^^,  the  yield  of  LBH  1400-1500  A 
/  =  2  fovYl,%,  the  yield  of  LBH  1650-1800  A 

i  =  3  for£o^,  the  characteristic  energy  of  the  incident  proton  spectrum 
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The  uncertainties  of  the  calculated  quantities  are  estimated  using  Equations  (3)  and  (4). 

In  addition,  READCO  reads  coefficients  and  parameters  for  ELAYER,  which  calculates 
N„E  and  hE  and  their  uncertainties.  This  model  is  described  in  Section  2.1 1. 

m  m 


2.4  SUBROUTINE  READGP 

SUBROUTINE  READGP  is  called  from  the  main  program  with  five  arguments; 

F10P7;  the  solar  activity  index  ,  [REAL,  Solar  Flux  units] 

AP:  the  magnetic  activity  index  [REAL,  unitless] 

KP:  the  magnetic  activity  index  [REAL,  unitless] 

QEUV:  the  solar  activity  index  determined  from  the  SSUSI  daytime  algorithm 
[REAL,  erg  cm s  '] 

LUN:  the  logical  unit  number  of  the  file  containing  the  first  three  variables  [INTEGER] 

READGP  reads  the  geophysical  parameters  from  a  file.  In  the  operational  version,  this  should  be 
replaced  with  a  call  to  the  appropriate  data  in  the  Environmental  Database  at  the  Space  Forecast 
Center.  However,  the  Q^^y  parameter  should  be  passed  in  some  manner  from  the  Daytime 
Algorithm.  In  the  current  version  of  this  algorithm,  these  parameters  are  used  only  to  estimate  the 
solar  contribution  to  the  E-layer  electron  density. 


2.5  SUBROUTINE  READIMG 

SUBROUTINE  READIMG  is  called  from  the  main  program  with  one  argument: 

LUN;  the  logical  unit  number  of  the  file  containing  the  images  [INTEGER] 

READIMG  reads  the  pixel  parameters  and  stores  them  in  COMMON  /IMAGE/  (Section  2.1). 

In  the  operational  version,  this  subroutine  must  be  modified  to  read  image  files  in  the  format 
produced  by  APL  algorithms. 


2.6  SUBROUTINE  GEOCRNA 

SUBROUTINE  GEOCRNA  is  called  from  the  main  program  with  5  arguments: 

A:  The  constant  term  in  the  functional  form  of  the  geocoronal  background 

[REAL,  Rayleighs]. 

B:  The  coefficient  of  the  cosine  solar  zenith  angle  term  in  the  functional  form  of 

the  geocoronal  background  [REAL,  Rayleighs]. 
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AVAR:  The  variance  of  A  [REAL,  Rayleighs^] 

BVAR:  The  variance  of  B  [REAL,  Rayleighs^] 

ABVAR:  the  covariance  of  A  and  B  [REAL,  Rayleighs^] 

GEOCRNA  also  has  access  to  the  contents  of  COMMON  /IMAGE/. 

GEOCRNA  does  not  actually  subtract  the  geocoronal  Lyman  a  background  from  the 
Lyman  a  image,  but  merely  returns  coefficients  that  allow  the  subtraction  to  be  done  later.  The 
algorithm  is  based  on  the  assumption  that  the  geocoronal  background  may  be  represented  by  the 
simple  functional  form 

I^  =  A  +  Bcosx  (6) 

where  4  is  the  background  intensity  in  Rayleighs  and  x  is  the  solar  zenith  angle.  The  algorithm 
assumes  that  the  auroral  boundaries  have  already  been  identified  and  that  each  pixel  is  flagged  so 
that  it  can  be  identified  as  either  auroral  or  non-auroral.  Only  non-auroral  pixels  are  used  to 
determine  A  and  B  by  means  of  a  straightforward  least  squares  fit.  The  fitting  process  produces 
estimates  of  the  variances  of  A  and  B,  as  well  as  the  covariance  of  the  two  with  each  other,  and 
these  are  used  to  estimate  the  “noise  level”  of  the  image.  These  are  returned  to  the  calling 
program  for  use  in  calculating  the  “noise  level,”  as 

=  2^0  ^  +a  \  cos^  X  +  2c  cos  x  (7) 

where  is  the  variance  of  A,  o^is  the  variance  of  B,  and  is  the  covariance  of  A  and  B. 
(Whenever  the  difference  between  the  measured  intensity  and  the  modeled  background  intensity  is 
less  than  the  noise  level,  the  difference  is  set  to  zero.) 


2.7  SUBROUTINE  DAYGLO 

Like  GEOCRNA,  SUBROUTINE  DAYGLO  determines  coefficients  that  allow  the 
subtraction  of  the  dayglow  from  the  two  LBH  images.  It  is  called  from  the  main  program  with 
ten  arguments. 

A1450:  the  constant  term  of  the  LBH  1400-1500  A  dayglow  expression  [REAL, 

Rayleighs] 

B1450:  the  coefficient  of  the  cosine  term  of  the  LBH  1400-1500  A  dayglow  [REAL, 

Rayleighs] 

VA1450:  the  variance  of  A 1450  [REAL,  Rayleighs^] 

VB1450:  the  variance  of  B1450  [REAL,  Rayleighs^] 

VC1450;  the  covariance  of  A1450  and  B1450  [REAL,  Rayleighs^] 

A1725;  the  constant  term  of  the  LBH  1650-1725  A  dayglow  expression  [REAL, 

Rayleighs] 
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B 1 725 :  the  coefficient  of  the  cosine  term  of  the  LBH  1 650- 1 725  A  dayglow  expression 

[REAL,  Rayleighs] 

VA1725:  the  variance  of  A1725  [REAL,  Rayleighs^] 

VB 1 725 :  the  variance  of  B 1 725  [REAL,  Rayleighs^] 

VC1725:  the  covariance  of  A1725  and  B1725  [REAL,  Rayleighs^] 


In  addition,  DAYGLO  has  access  to  the  contents  of  /IMAGE/. 


This  algorithm  is  similar  to  GEOCRNA  except  that  instead  of  fitting  all  non-auroral 
pixels,  it  fits  only  day  pixels.  The  fitting  function  similar  to  the  one  used  in  GEOCRNA; 


Ch(x) 


(8) 


where  X  is  either  145  nm  for  LBHl  or  1 72.5  nm  for  LBH2,  %  is  the  solar  zenith  angle,  and  Ch(x) 
is  the  Chapman  grazing  incidence  function,  which  we  approximate  by  [Rishbeth  and  Garriott, 
1969] 


Ch(x)s 


^^Tcxsinx 

l-erf(^ 

^jxcos'  x) 

^yTtxsinx 

1  +  erf(^ 

^Ixcos^  x) 

(9) 


wherex  =  -t-z)  /  .  We  treat  x  as  a  constant  with  a  value  of  720  (obtained  by  setting  z  = 

1 10  km  and  =  9  km).  We  approximate  the  error  function  as  [Abramowitz  and  Stegun,  1964]: 


with 


and 


erf(w)  =  \-{a^t  +  a.^t'  -\-a^t^  +aj'*  +a^t^^e~'‘^  +8(m)  (10) 


1 

r  = - 

\  +  pu 

/?  =  0.3275911 

a,  =  0.254829592 

a2= -0.284496736 

a  =  1.421413741 

a,  =  -1.453152027 

a^=  1.061405429 

(For  X  less  than  about  60°,  Ch(x)  is  well  approximated  bysecx  •) 


The  variance  of  is 

cos-x+2a^^B^cosx  (11) 

where  the  variances  and  covariances  of  and  are  determined  by  the  least  squares  fitting 
process. 
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2.8  SUBROUTINE  AURORA 


SUBROUTINE  AURORA  is  the  primary  subprogram  of  the  auroral  algorithm,  it  loops 
over  pixels,  calling  SUBROUTINE  PROTON,  SUBROUTINE  ELECTRN,  and  SUBROUTINE 
ELAYER  as  needed.  It  is  called  from  the  main  program  with  27  arguments. 


A: 

B: 

AVAR: 

BVAR: 

ABVAR: 

A1450: 

B1450: 

VA1450: 

VB1450: 

VC1450: 

A1725; 

B1725: 

VA1725: 
VB1725: 
VC  1725: 
EOE: 

QOE: 

EOP: 

QOP. 

NME: 

HME: 


UQOE: 

UEOP: 

UQOP: 


the  constant  term  of  the  geocorona  background  [REAL,  Rayleighs] 

the  coefficient  of  the  cosine  solar  zenith  angle  term  of  the  geocorona 

background  [REAL,  Rayleighs] 

the  variance  of  A  [REAL,  Rayleighs^] 

the  variance  of  B  [REAL,  Rayleighs^] 

the  covariance  of  A  and  B  [REAL,  Rayleighs^] 

the  constant  term  of  the  LBH  1400-1500  A  dayglow  expression  [REAL, 

Rayleighs] 

the  coefficient  of  the  cosine  term  of  the  LBH  1400-1500  A  dayglow  expression 
[REAL,  Rayleighs] 

the  variance  of  A1450  [REAL,  Rayleighs^] 

the  variance  of  B 1450  [REAL,  Rayleighs^] 

the  covariance  of  A1450  and  B1450  [REAL  ARRAY,  Rayleighs^] 

the  constant  term  of  the  LBH  1650-1725  A  dayglow  expression  [REAL, 

Rayleighs] 

the  coefficient  of  the  cosine  term  of  the  LBH  1650-1725  A  dayglow  expression 
[REAL,  Rayleighs] 

the  variance  of  A 1725  [REAL,  Rayleighs^] 

the  variance  of  B 1725  [REAL  ARRAY,  Rayleighs^] 

the  covariance  of  A1725  and  B1725  [REAL  ARRAY,  Rayleighs^] 

the  array  containing  the  deduced  electron  characteristic  energy  [REAL 

ARRAY,  keV] 

the  array  containing  the  deduced  electron  energy  flux  [REAL  ARRAY, 
erg  cm s  '] 

the  array  containing  the  deduced  proton  characteristic  energy  [REAL  ARRAY, 
keV] 

the  array  containing  the  deduced  proton  energy  flux  [REAL  ARRAY, 
erg  cm ^  s  '] 

the  array  containing  the  deduced  peak  E-layer  density  [REAL  ARRAY,  cm'^] 
the  array  containing  the  deduced  E-layer  height  [REAL  ARRAY,  km] 
the  array  containing  the  uncertainty  in  the  deduced  electron  characteristic 
energy  [REAL  ARRAY,  keV] 

the  array  containing  the  uncertainty  in  the  deduced  electron  energy  flux  [REAL 
ARRAY,  erg  cm s  '] 

the  array  containing  the  uncertainty  in  the  deduced  proton  characteristic  energy 
[REAL  ARRAY,  keV] 

the  array  containing  the  uncertainty  in  the  deduced  proton  energy  flux  [REAL 
ARRAY,  erg  cm ^  s  '] 
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UNME: 


the  array  containing  the  uncertainty  in  the  deduced  peak  E-layer  density 
[REAL  ARRAY,  cm^] 

UHME:  the  array  containing  the  uncertainty  in  the  deduced  E-layer  height  (REAL 

ARRAY,  km] 


AURORA  also  has  access  to  the  contents  of  /IMAGE/. 

AURORA  loops  over  all  pixels.  If  a  pixel  is  tagged  as  auroral,  then  it 

(1)  subtracts  the  geocoronal  and  dayglow  backgrounds  and  estimates  the  uncertainty  of 
the  result, 

(2)  calls  PROTON  with  MODE=0, 

(3)  calls  ELECTRN,  and 

(4)  if  and  only  if  QOE  is  small  enough  (QOE  <  0.001  QOP),  it  calls  PROTON  with 
(MODE-1). 


Otherwise,  it  merely  skips  the  pixel. 

The  background  subtraction  and  uncertainty  estimation  proceeds  as  follows.  For  the 
geocorona 


*  geocorona 


=  A  +  Bcosx 


(12) 


and 


geocorona 


=  cl+alcos-x+2o^^BCOsx 


(13) 


The  proton  contribution  to  the  Lyman  a  intensity  is  estimated  as 


jip)  _  [■^1216  "  ^ geocorona 
•‘1216  ~ 


0, 


1216  ^  ^  geocorona  geocorona 

Otherwise 


(14) 


where  7,2,6  is  the  measured  intensity  at  121.6  nm. 
The  uncertainty  in  is  estimated  by 


^  1216,/j  ^  ]2\6  ^  geocorona 


(15) 


For  the  dayglow 


D 

"  "  Ch(x) 


(16) 
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where  X  is  either  145  nm  for  LBHl  or  172.5  nm  for  LBH2,  %  is  the  solar  zenith  angle,  and  Ch(x) 
is  the  Chapman  grazing  incidence  function  described  in  Section  2.5  above. 


The  variance  of  is 

cos'  X  +  2o cos X 


(17) 


As  for  the  geocoronal  subtraction,  the  auroral  (proton  +  electron)  contribution  tol^,  the  measured 
intensity  at  wavelengthX,  is  estimated  by 


4'"’ 


0,  otherwise 


The  uncertainty  in  is  estimated  from 


a 


2 

X,e+  p 


(18) 


(19) 


2.9  SUBROUTINE  PROTON 

SUBROUTINE  PROTON  is  called  from  AURORA  with  15  arguments: 

11216;  the  Lyman  a  intensity  after  background  subtraction  [REAL,  Rayleighs] 

11450:  the  LBH  1  (1400  -  1500  A)  intensities  after  dayglow  subtraction  [REAL, 

Rayleighs] 

11725:  the  LBH  2  (1650  -  1800  A)  intensities  after  dayglow  subtraction  [REAL, 

Rayleighs] 

U1216:  the  uncertainty  in  the  Lyman  a  intensity  after  background  subtraction  [REAL, 

Rayleighs] 

U1450:  the  uncertainty  in  the  LBH  1  intensities  after  dayglow  subtraction  [REAL, 

Rayleighs] 

U1725;  the  uncertainty  in  the  LBH  2  intensities  after  dayglow  subtraction  [REAL, 
Rayleighs] 

EOE;  the  characteristic  energy  of  the  precipitation  electrons.  [REAL,  keV] 

NOTE;  If  MODE  (see  below)  is  0,  then  this  parameter  is  unused. 

QOE:  the  energy  flux  of  precipitating  electrons.  [REAL,  erg  cm'^  s  '] 

NOTE:  If  MODE  is  0,  then  this  parameter  is  unused. 

UEOE;  the  uncertainty  in  EOE  [REAL,  keV] 

UQOE:  the  uncertainty  in  QOE  [REAL,  erg  cm'^  s  '] 

MODE:  an  integer  variable  indicating  the  "mode"  in  which  PROTON  is  to  calculate  the 
proton  precipitation  parameters.  [INTEGER] 

EOP:  the  characteristic  energy  of  the  precipitation  protons.  [REAL,  keV] 

NOTE;  If  MODE  is  0,  then  this  parameter  is  set  to  8  keV. 

QOP:  the  energy  flux  of  precipitating  protons.  [REAL,  erg  cm  '^  s  '  ] 
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UEOP:  the  uncertainty  in  EOP  [REAL,  keV] 

UQOP:  the  uncertainty  in  QOP  [REAL,  erg  cm'^  s  '] 

PROTON  operates  on  a  single  pixel  all  of  the  arguments  are  scalars.  There  are  two  possible 
modes;  MODE  =  0  when  the  electron  contribution  is  unknown,  and  MODE  =  1  when  the 
electron  contribution  is  known. 


For  MODE  =  0,  PROTON  estimates  the  incident  proton  energy  flux,  by  assuming  that 
£0  p  =  8  keV.  The  uncertainty  in  ^  is  set  to  4  keV  (50%  of  the  assumed£o  p). 

For  MODE  =  1,  PROTON  subtracts  the  electron  contribution  to  the  two  LBH  bands 
using  EOE  and  QOE.  It  then  estimates  both  ^  from  the  ratio  of  the  two  LBH  bands  and 

from  the  Lyman  a  intensity. 

As  a  practical  matter,  there  is  no  way  to  estimate  the  electron  contribution  independent  of 
the  proton  contribution  using  SSUSI  data  alone.  Therefore,  MODE  =  1  is  used  only  when  the 
electron  contribution  is  believed  to  be  zero  or  negligible.  The  normal  calling  sequence  is 

1 .  Call  PROTON  with  MODE  =  0  to  get  an  estimate  of  the  proton  flux. 

2.  Call  ELECTRN  to  get  an  estimate  of  the  electron  flux. 

3 .  If  and  only  if  the  estimate  of  the  electron  energy  flux  is  less  than  1 0'^  times  the  proton 
flux.  Call  PROTON  with  MODE  =  1  and  QOE  =  0.0  to  get  a  refined  estimate  of  the 
proton  energy  flux  and  mean  energy. 

The  proton  energy  flux  is  estimated  from  an  expression  of  the  form 

lip) 

(20) 

^1216 


where  is  the  Lyman  a  yield  of  a  1  erg  cm  '^  s  '  proton  aurora,  and 


'1216 


Ageocorona) 
■'1216  > 


otherwi  se 


(21) 


The  yield  is  estimated  from  the  series 


When  MODE  =  0,  Ea  j,\s  set  to  8  keV.  Otherwise  (MODE  =  1),  E^^  ^is  determined  from  the  ratio 
of  the  1400-1500  A  LBH  and  1650-1800  A  LBH  intensities  through  an  expression  of  the  form 


15 


(23) 


(;’2)  ,  _(p2)  n 
0  ■'■‘'1  -^LSW 


n=0 


where 


and 


^LBH  ~ 


A''  =  4 


j{dayglow)  _  j(e) 


(24) 


(25) 


with  X  either  1450  A  or  1725  A.  (The  background  subtraction  is  assumed  to  have  been  carried 
out  before  the  call  to  PROTON) 


In  MODE  =  1,  the  electron  contribution  to  the  two  LBH  intensities  is  estimated  from 
expressions  of  the  form 


where 

n=0 

n^O 

Of  course,  in  normal  usage  will  be  zero  and  both  and  will  vanish.  As  delivered,  the 
FORTRAN  code  tests  QOE,  and  if  it  is  exactly  zero,  the  electron  contribution  is  set  to  zero 
without  evaluating  the  above  series. 


2./0  SUBROUTINE  ELECTRN 

SUBROUTINE  ELECTRN  is  called  from  AURORA  with  10  arguments; 

11450:  the  summed  intensities  in  the  LBH  1  wavelength  range  (1400-1500  A)  [REAL, 

Rayleighs] 

11725:  the  summed  intensities  in  the  LBH  2  wavelength  range  (1650-1800  A)  [REAL, 

Rayleighs] 

C1450;  the  actual  number  of  counts  due  to  the  LBH  1  intensity  [INTEGER,  unitless] 
Cl  725:  the  actual  number  of  counts  due  to  the  LBH  2  intensity  [INTEGER,  unitless] 
EOP:  the  characteristic  energy  of  the  precipitation  protons.  [REAL,  keV] 
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QOP:  the  energy  flux  of  precipitating  protons.  [REAL,  erg  cm'^  s'] 

EOE:  the  characteristic  energy  of  the  precipitation  electrons.  [REAL,  keV] 

QOE:  the  energy  flux  of  precipitating  electrons.  [REAL,  erg  cm'^  s  '] 

UEOE:  the  uncertainty  in  the  characteristic  energy  of  the  precipitation  electrons. 

[REAL,  keV] 

UQOE:  the  uncertainty  in  the  energy  flux  of  precipitating  electrons.  [REAL, 

erg  cm'^  s  '] 


Like  PROTON,  ELECTRN  operates  on  a  single  pixel,  so  all  of  the  arguments  are  scalars. 

ELECTRN  calculates  the  energy  flux  and  characteristic  energy  from  the  two  LBH  bands 
after  subtracting  the  proton  contribution.  If  the  proton  energy  flux  is  non-zero,  ELECTRN 

calls  the  FUNCTION  subprogram  PSEVAL  to  estimate  the  proton  contribution: 

with 

=  (28) 

rt=0 

=  (29) 

rt=0 


If  the  estimated  proton  contribution  is  greater  than  or  equal  to  the  observed  intensity  of  either 
LBH  feature,  the  electron  contribution  is  assumed  to  vanish.  Otherwise,  the  electron  contribution 
is  just  the  difference  between  the  observed  intensity  and  the  estimated  proton  contribution: 


0 


if  >  / 

^1450  —  ■'1450 


/<*) 

1/  if/'’^  </ 

L''1450  ■'1450  “  '1450  ^  '1450 


(30) 


/(«)  _ 
'l725  ~  ' 


0 


if/</’>  >  / 

“  'l725  -  'l725 


/  if/''^  <7 

'1450  '1450  “  '1725  ^  'l725 


(31) 


If  the  estimated  electron  contribution  to  either  LBH  feature  is  zero,  then  the  electron  energy  flux 
is  also  set  to  zero,  and  the  characteristic  energy  is  set  to  a  minimum  value  (500  eV). 

The  calculation  of  electron  precipitation  parameters  uses  the  following  relationship 
between  characteristic  energy,  and  the  ratio  of  the  two  intensities7|45o  and/j^j: 
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n=0 

(32) 

where 

/(e) 

tL‘)  _  ^1450 
-  j{e) 

^1725 

(33) 

Once  £0,  has  been  determined,  the  energy  flux  is  estimated  in  the  following 

manner; 

n=0 

(34) 

>1=0 

(35) 

give  the  yields  of  the  two  LBH  bands  for  an  electron  energy  flux  ofl  erg  cm' 
observed  intensity  to  the  yield  gives  the  energy  flux: 

^  s'' .  The  ratio  of  the 

/(«) 

WI450)  _  41450 

4  1450 

(36) 

/(«) 

Wm5)  _  4,725 
y(e) 

4 1725 

(37) 

Since  it  is  likely  that  and  will  differ,  ELECTRN  returns  the  value  determined  from 

the  feature  with  the  largest  comUitig  rate.  That  is 

ifC  >C 

^  “  '-^1450  — '^1725 

1/^(1725)  -r.^  ^ 

>  **  ^^1450  ^ '^1725 

(38) 

so  that  is  determined  from  the  feature  that  has  the  lowest  uncertainty  due  to  counting 
statistics. 
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2.11  SUBROUTINE  ELAYER 


SUBROUTINE  ELAYER  is  called  from  the  main  program  with  1 5  arguments: 


SZA: 

EOE: 

QOE: 

EOF; 

QOP: 

QEUV: 

UEOE; 

UQOE: 

UEOP: 

UQOP: 

UQEUV: 

NME: 

HME: 

UNME: 

UHME: 


solar  zenith  angle  [REAL,  degrees] 

the  characteristic  energy  of  the  precipitating  electrons  [REAL,  keV] 
the  energy  flux  of  the  precipitating  electrons  [REAL,  erg  cm'^  s'] 
the  characteristic  energy  of  the  precipitating  protons  [REAL,  keV] 
the  energy  flux  of  the  precipitating  protons  [REAL,  erg  cm  s'] 
the  solar  energy  flux  determined  by  the  daytime  algorithm  [REAL, 
erg  cm'^  s''] 

the  uncertainty  in  the  characteristic  energy  of  the  precipitating  electrons 
[REAL,  keV] 

the  uncertainty  in  the  energy  flux  of  the  precipitating  electrons  [REAL, 
erg  cm'^  s  '] 

the  uncertainty  in  the  characteristic  energies  of  the  precipitating  protons 
[REAL,  keV] 

the  uncertainty  in  the  energy  fluxes  of  the  precipitating  protons  [REAL, 
erg  cm  s''  ] 

the  uncertainty  in  the  solar  activity  index  determined  by  the  daytime 
algorithm  [REAL,  erg  cm'*  s  '] 

the  electron  density  at  the  peak  of  the  E-layer  [REAL,  cm  '^] 
the  height  of  the  peak  of  the  E-layer  [REAL,  km] 
the  uncertainty  in  the  peak  electron  density  [REAL,  cm'^] 
the  uncertainty  in  the  height  of  the  peak  [REAL,  km] 


Like  PROTON  and  ELECTRN,  ELAYER  operates  on  one  pixel  at  a  time,  so  all  of  the  arguments 
are  scalars. 


The  first  step  is  the  determination  of  the  height  of  the  peak  ionization  due  to  the  two 
auroral  ionization  sources  (electrons  and  protons): 


“ICy 


b  =  e,p 


(39) 


where  =  1  keV  and  =  4  keV.  The  second  step  is  the  determination  of  the  peak 
production  rate  for  the  auroral  ionization  sources  for  an  incident  flux  of  1  erg  cm'^  s  ' : 


f 

•og,o  ■ 
V 


=  Ir 


«>) 


log 


( E  ^ 

^O.b 
>0  rib) 

J 


h  =  e,p 


(40) 
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The  next  step  is  the  representation  of  all  three  production  rate  profiles  by  Chapman  layers. 
We  first  define  a  generalized  Chapman  function  S{u,\]/)  as 

=  exp[l-i/-e'"Ch(v|/)]  (41) 

where  Ch(v|/)  is  the  Chapman  grazing  incidence  function  described  in  Section  2.5.  We  write 


a 


erg  cm s  ' 


b  =  e,p 


(42) 


for  protons  and  electrons,  and 


a 


1  erg  cm  s 


-pT'''\0)S{ueuy^x) 


(43) 


for  photoionization.  In  the  above  expressions 


b  =  e,p 


and 


z-K 


(44) 


(45) 


For  protons  and  electrons,  /r^*Mepends  on  the  characteristic  energy,  and  is  given  by 
Equations  (33)  and  (32),  respectively.  For  photoionization,  //^is  just  the  height  of  the  peak  E- 
layer  ionization  at  the  subsolar  point  (~108  km),  and  is  just  the  neutral  scale  height  in  the  E 
region  (~9  km).  For  protons  and  electrons,  determining  requires  a  little  more  effort. 

First  we  note  that  although  the  column  ionization,  Q]'’\  is  independent  ofE^^,  the  peak 
ionization,  pl^\  is  not.  This  implies  that  the  shape  of  the  ionization  profile  changes  with 
characteristic  energy.  In  the  Chapman  function,  the  shape  is  controlled  by  the  scale  height, 
which  may  be  determined  in  the  following  manner. 


h,qj:' 

1  erg  cm s  ' 


J  S{u^,0)du, 

¥—tt> 


(46) 


Since 


f  S{u,\\i)du  =  e{  e  W  = 

J-QO  Jo 


e 

Ch(\|/) 


(47) 
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(where  e  is  the  base  of  the  natural  logarithms,  not  the  charge  on  an  electron),  we  have 


(b) 


1  erg  cm s  ' 


(48) 


or 


1(6) 


_  1  erg  cm  s  Q] 

“  a 


(49) 


For  electrons  and  protons  the  total  column  ionization  rate  Q\’’^  in  cm'^  s  '  is  related  to  the 
incident  energy  flux  of  that  species  Q^,  in  erg  cm'^  s'  by 


a*'' = p.a 


(50) 


The  proportionality  constant  p*  includes  the  conversion  from  erg  to  eV  and  the  mean  “eV  per  ion 
pair”: 


_ 1 _ 

43.73  ev/ ion  pair  X  1.6022  xlO'''erg/eV 
_ 1 _ 

27  ev/ion  pair  x  1.6022  x  10''^ erg/ eV 


w  1.427x10'°  ion  pairs /erg 
«  2.3  X 10'°  ion  pairs /erg 


The  “eV  per  ion  pair”  for  electrons  differs  from  the  traditional  “35  eV  per  ion  pair”  because  a 
significant  fraction  of  the  incident  electron  flux  is  backscattered  out  of  the  atmosphere  before  it 
can  do  any  appreciable  ionization.  H^\s  be  determined  from 


fib{Eo,)  = 


Pt  x(l  erg  cm"  s  ') 


(51) 


The  altitude  distribution  of  ionization  is  now  completely  specified  for  all  three  ionization  sources. 

The  composite  electron  density  profile  is  calculated  using  an  effective  recombination  rate, 
a{z),  which  is  a  function  of  altitude: 


,,(.)=  jSEZMIS  (52) 

where  «^(z)  is  the  electron  density.  By  squaring  both  sides  of  Equation  (52)  and  taking  a 
derivative  with  respect  to  z,  one  may  obtain  an  condition  for  the  maximum  electron  density: 
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(53) 


or  more  simply 

1  dpf  _  1  da 
dz  a  dz 

where 


A(^)  = 

b 

A  somewhat  crude  but  serviceable  model  for  a{z)  is 


a(")  =  1 


a 


0’ 


r 


ttoCxp 


V 


z<z„ 


z>z„ 


(54) 


(55) 


(56) 


Consequently,  if  z  <z^,  h^E  is  the  same  as  the  height,  /?„,  of  the  peak  in  the  total  ionization,  p,  . 
In  general,  h^Eis  determined  by 


Pt 


Pi 


h„<z„ 


(57) 


Unfortunately,  it  is  not  possible  to  derive  an  expression  for  h^E  in  closed  form  unless  a  single 
ionization  source  dominates.  In  that  case 


1  dp'^  1  15%,0 

(t)  '  =  ^  = - E —  ionization  source) 

p]  ’  dz  H,S  dii,  H, 


(58) 


From  Equation  (57) 


f  *1*’. 


hE=< 


1  + 


a  J 


(single  ionization  source)  (59) 


When  two  or  more  of  the  ionization  sources  produce  comparable  ionization  rates,  h^E 
must  be  determined  using  some  sort  of  numerical  root  finding  algorithm.  We  have  decided  to  use 
the  simplest  such  algorithm:  simply  tabulate  n^{z)  between  90  and  160  km  at  Az  =  5  km  intervals 
and  find  the  largest  “internal”  maximum.  (Specifically,  we  reject  a  maximum  at  the  ends  of  the 
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interval.)  If  there  is  no  internal  maximum,  then  we  set  h^E  to  1 10  km.  Although  this  method  is 
rather  crude,  it  has  the  virtue  of  requiring  only  14  evaluations  of  w,(r).  Given  the  unpredictable 
and  variable  nature  of  the  auroral  EDP  in  the  £'-layer,  there  is  no  guarantee  that  a  more 
sophisticated  root  finding  algorithm  would  converge  in  a  reasonable  number  of  iterations.  Using 

Az 

the  simple  method,  the  uncertainty  in  is  no  more  than  ±2.5  km  (“~)  unless  no  internal 

maximum  is  found.  In  that  case,  however,  the  auroral  contribution  to  the  E-layer  must  be 
unimportant,  so  the  uncertainty  in  h^E  is  also  unimportant. 

Once  hE  has  been  determined,  N„E  is  calculated  from 

m  ^  m 


(60) 


where  p{z)  =  ^^p\^\z)  is  the  total  ionization  rate  at  altitude  z. 

b 

Because  of  the  inherent  nonlinearity  of  these  formulae,  error  analysis  for  N^E  and  h^E  is 
quite  complicated.  As  a  practical  matter,  we  assume  that  the  uncertainty  in  h^E  calculated  from 
Equation  (57)  is  no  more  than  the  error  in  h„,E  determined  by  the  tabulation  method  above.  Then 


a  conservative  estimate  of  the  uncertainty  in  h^,E  is  v2—  =  3.5  km,  and  the  uncertainty  in 


N„,E  =  nXKE)  is  simply 


1^1’  I 

N„E 


(61) 


where  we  have  assumed  that  the  errors  in  p  and  a  uncorrelated.  However,  this  does  not  take  into 
account  the  error  in  N„E  due  to  errors  in  h„E.  These  are  hard  to  estimate  in  this  calculation 

m  m 

method,  but  we  will  assume  that  they  are  no  larger  than  the  uncertainty  estimate  of  Equation  (61). 
In  that  case,  a  conservative  estimate  of  the  uncertainty  in  N^E  is 


^  N„E  hp  Og 

N„,E  2  ip^ 

The  uncertainty  in  p  may  be  estimated  from 


(62) 


(63) 
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where  the  uncertainties  in  are  to  be  evaluated  at  z  =  h^E.  These  uncertainties  depend  on  the 
uncertainty  in  the  Chapman  function: 


with 


[e-'‘Ch{x)-\fol 


(64) 


(65) 


Similarly,  the  uncertainty  in  a  may  be  estimated  from 


(66) 


where  the  a  parameters  are  assumed  to  be  uncorrelated. 

SSUSIAE  evaluates  uncertainties  in  each  quantity  as  it  is  calculated  and  passes  the 
estimated  uncertainty  to  the  next  subprogram  that  needs  it.  For  mathematical  completeness,  an 
error  analysis  giving  the  uncertainties  in  N^E  and  h^E  in  terms  of  each  primitive  parameter  is 
provided  by  the  Appendix. 
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2.12  SUBROUTINE  WIMAGE 

SUBROUTINE  WIMAGE  is  called  from  the  main  program  with  14  arguments; 

NPIXEL:  The  number  of  pixels  in  the  image  [INTEGER] 

EOE:  the  array  that  contains  the  characteristic  energies  of  the  precipitating  electrons 

as  deduced  from  the  image  intensities  [REAL  ARRAY(NPIXEL),  keV] 

QOE:  the  array  that  contains  the  energy  fluxes  of  the  precipitation  electrons  as 

deduced  from  the  image  intensities  [REAL  ARRAY(NPIXEL),  erg  cm'^  s  '  ] 
EOF:  the  array  that  contains  the  characteristic  energies  of  the  precipitating  protons  as 

deduced  from  the  image  intensities  [REAL  ARRAY(NPIXEL),  keV] 

QOP;  the  array  that  contains  the  energy  fluxes  of  the  precipitation  protons  as 

deduced  from  the  image  intensities  [REAL  ARRAY(NPIXEL),  erg  cm'^  s  '  ] 
NME;  the  array  that  contains  the  peak  electron  density  of  the  E  layer  (N^E)  as 
deduced  from  EOE,  QOE,  EOF,  and  QOE  [REAL  ARRAY,  cm’'] 

HME:  the  array  that  contains  the  height  of  the  E-layer  [REAL  ARRAY,  km] 

UEOE;  the  array  that  contains  the  uncertainties  in  the  characteristic  energies  of  the 

precipitating  electrons  [REAL  ARRAY(NFIXEL),  keV] 

UQOE:  the  array  that  contains  the  uncertainties  in  the  energy  fluxes  of  the  precipitation 

electrons  [REAL  ARRAY(NFIXEL),  erg  cm'^  s  '  ] 

UEOF:  the  array  that  contains  the  uncertainties  in  the  characteristic  energies  of  the 

precipitating  protons  [REAL  ARRAY(NFIXEL),  keV] 

UQOF:  the  array  that  contains  the  uncertainties  in  the  energy  fluxes  of  the  precipitation 

protons  [REAL  ARRAY(NFIXEL),  erg  cm'^  s  ’  ] 

UNME:  the  array  that  contains  the  uncertainties  in  the  peak  electron  density  of  the  E 

layer  (N„,E)  [REAL  ARRAY,  cm”] 

UHME:  the  array  that  contains  the  uncertainties  in  the  height  of  the  E-layer  [REAL 

ARRAY,  km] 

LUN:  the  logical  unit  number  of  the  file  where  the  "image"  is  to  be  written. 

WIMAGE  writes  an  "image"  file  in  the  same  format  as  the  original  image  file  except  that 
precipitation  parameters  (E^^,  Q^,  E^^,  Q^,  N^E,  and  h^E)  and  their  estimated  uncertainties, 
rather  than  intensities,  are  being  stored. 

WIMAGE  is  provided  solely  to  allow  AURORA  to  operate  as  a  stand-alone  algorithm.  We 
assume  that  the  final  package  of  SSUSI  algorithms  will  have  its  own  I/O  module(s)  and  will 
dispense  with  WIMAGE. 
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APPENDIX;  h^E  and  N^E  error  propagation  analysis  in  terms  of  primitive  variables 

Both  h„E  and  N„E  depend  on  six  “arguments”:  Q^,  Q^,  Q^y,  E^  „  E^  ^,  and  x  •  (In  the 
following  error  propagation  analysis,  the  uncertainty  in  x  is  assumed  to  be  negligible.)  In 
addition,  both  depend  on  a  number  of  other  parameters  determined  by  fitting  the  results  of 
computer  modeling  with  analytic  functions.  These  parameters  are 


Parameter _ Description _ Equation(s) 


Pi‘’\  n=0,l,2,3;  b=e,p 

coefficients  in 

(39) 

K 

height  of  peak  photoionization  at  x  =  0 

(43) 

neutral  scale  height  at 

(43) 

w=0,l,2;  b=e,p 

coefficients  in  p^^^{Eq^) 

(40) 

p'^{^::}),b=e,p 

peak  impact  ionization  rate  at  the  reference 
energy  of  (sometimes  denoted  by  p^^l^f) 

(40) 

pT\o) 

peak  photoionization  at  x  =  0  (sometimes 
denoted 

(43) 

Pft,  b=e,p 

proportionality  constant  between  andQ^ 

after  (50) 

the  altitude  at  which  a(z)  changes  from  constant 
to  and  exponential  function 

(56) 

Cto 

a(^a) 

(56) 

Ha 

the  scale  height  of  a{z)  for  z>z^ 

(56) 

The  general  equation  for  calculating  the  uncertainties  of  h^E  and  N^^E  are  derived  here. 
Explicit  expressions  in  terms  of  the  above  parameters  follow  the  general  equations. 

Both  quantities  may  be  thought  of  as  functions  of  the  form 


f{^Qe->QpyQEVVi^0,ei^0.p'>  X>*) 


where  a  is  a  vector  whose  elements  are  all  the  parameters  used  in  the  E  layer  model,  including 
P!.'’’.  y!,*’.  Pe>  Pp.  and  The  uncertainty  in /may  be 

estimated  as  [Bevitigton  and  Robinson,  1992] 


^  f-^Q, 


^  df^ 
\^QeJ 


+a; 


^  df^ 


+a 


Qkut 


df 


dQ, 


+a 


V  '^iCEuy , 


1  f  5/  Y  2 


y^^o.p  j 


^  df  Y 


(Al) 
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where  we  have  assumed  that  the  uncertainty  in  %  is  negligible  and  that  the  six  arguments  are 
uncorrelated  with  each  other  and  with  the  parameters  a.  Of  course,  many  of  the  parameters  are 
also  uncorrelated  with  other  parameters,  so  the  covariance  matrix  a  has  many  zero  elements. 
It  must  also  be  remembered  that  depends  explicitly  on  h^E,  so  the  partial  derivatives  with 
respect  to  the  arguments  and  parameters  must  include  the  implicit  dependences  represented  by 
KE. 


If  h^E  were  determined  numerically  from  Equation  (57),  the  partial  derivatives  would 
have  to  be  calculated  implicitly.  Let 


g(Qe,  Qp^  QeUV^  Eo.e,  ^ ^ 

so  that  h^E  is  determined  implicitly  by 

s{Q,,  QpyQEUEyEo.ey  Eo.pyXyK,Ey»)  =  0 


(A2) 


(A3) 


The  partial  derivatives  required  for  the  error  propagation  analysis  are  obtained  as  follows. 

dg  _dg  ^  dg  dh„,E  _  ^ 
dh,  db,  dK,E  db, 


(A4) 


or 


\^KE)  db, 


(AS) 


where  b^.  can  be  any  of  the  arguments  (O,  E,  etc.)  as  well  as  any  parameters  represented  by  a. 
Taking  into  account  the  definition  of  g,  we  may  obtain 


dh„E 

db, 


dbjdz 


_L^ 

H^db, 


iz=h„E 


^  1  Sp. 

“ 

dz^ 

H.  Sz 

z=h,E_ 

(A6) 


which  applies  to  all  arguments  and  to  all  parameters  used  in  calculating  p.(z)  -  that  is  all 
parameters  except  a^, ,  ,  and  .  Explicit  expressions  for  the  partial  derivatives  in  Equation 

(A6)  are  given  below,  as  are  the  partial  derivatives  with  respect  to  the  three  a(z)  parameters. 
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The  derivatives  required  for  the  error  propagation  analysis  are  more  straightforward; 


dN„E 

db, 


1 

m 


Piifi.E) 


db, 


1  5p, 

aa 

^KE 

z-KE  ^  {KE)  ^ 

7=/l.£'y 

db, 

a\h„E) 


(A7) 


or 


^EfmE  jq  £ 

15/)^  1  5a  ^ 

1  _L 

'  \  dp^  1  5a^ 

dKE 

m 

56* 

56*  a56*J 

1  T 

'z=h,E 

Pi  dz  a  dz  J 

J 

However,  from  Equation  (54)  we  see  that  the  second  term  in  the  square  brackets  vanishes  at  z  = 
h^E  so  the  final  result  is 


SN„,E 

db, 


1  dPi 
.Pi 


1  5a  ^ 


(A8) 


The  partial  derivatives  needed  for  equation  (A8)  are  given  below.  N^E  is  also  an  explicit 
function  of  a(z)  and  its  parameters,  ag,  z„,  and  These  partial  derivatives  are  also  given 
below. 


The  explicit  expression  for  the  uncertainty  defined  by  Equation  (Al),  taking  into  account 
uncorrelated  parameters,  is: 


^  f~^Q. 


-fa 


Qp 


-fo 


Qew 


df 


5Q 


+a: 


\^\Ceuv  J 


+a1 


df 


3  3 


+ZZ^'p(')p(r 


^  df 


n=0  r/=0 


v3P'.'’y 


JL 


3  3 


2  2 


(«) 

j 

/  ar  ^  2  2 


+ZZ'’ 


n=0  tf-0 


JL 

3p': 


(.p) 


n=0  rf=0 


(e) 


df 


/.=0n’=0 


_  ,  -  1  H  t  if 


'  df 


(/») 


JL 


+o 


AK. 


+  CJ‘ 


^  df^ 


df 


+o 


n'*) 

Pm/tf 


df 


+o 


PmjtJ 


df 


(A9) 


J 
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J.Buy) 

Pm.O 


df 


dp] 


(Euy) 


+a 


p. 


K. 

dpe 


+ct; 


JL 

Jp.y 


JL 


+o; 


v^oy 


+a1 


JL 


a  / 
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Explicit  expressions  for  each  of  these  derivatives  are  given  below.  In  most  cases,  the  parameters 
appear  in  only  a  few  terms,  so  that  the  individual  partial  derivatives  are  relatively  simple.  Since  it 
is  essential  to  the  hJE  error  propagation  analysis,  the  partial  derivative  of  will  be 

evaluated  first. 


dhE  dz 

tn 


dz 


\z=h,E 


(AlO) 


where,  as  usual,  =  0  for  =  e  or  p,  and  =  % .  Evaluating  the  derivative  of  S  yields 


|-S(!/,v|/)  =  [«^‘Ch(i);)-l]S(K,v)  (All) 

which  appears  in  most  of  the  partial  derivatives  that  must  be  evaluated.  Also  important  is  the 
partial  derivative  of  dp^  /  dz : 


d 

dp, 

dhE 

tn 

dz 

dz^ 

dz 


2=h^E 


which  requires  the  evaluation  of  d^ SI diP' : 

Ch^(vi/ )  -  3^'"  Ch(v)/ )  +  l^z/.vi/ ) 

dir  *■ 


(A12) 


(A13) 


dh  E  dN  E 

The  general  expressions  for  — —  and  — —  are  given  by  Equations  (A6)  and  (A8)  above.  The 

d\  db^ 

expressions  are  in  terms  of  partial  derivatives  of  /?,  with  respect  to  b^,  which  may  represent  either 
an  argument  or  a  parameter.  It  is  these  partial  derivatives  that  are  evaluated  below. 

NOTE:  In  both  Equations  (A6)  and  (A8),  the  dependence  on  is  taken  into  account 
explicitly,  so  in  deriving  the  following  expressions  p^  was  treated  as  a  function  of  z  while  the 
argument  and  parameter  derivatives  were  taken,  and  then  the  result  was  evaluated  at  2  =  h^E. 


We  begin  by  evaluating  the  partial  derivatives  of  /?,  and  dp^  /  dz  with  respect  to  the  five 
relevant  arguments  (errors  in  %  being  ignored): 


ft: 


^  =  b  =  e.p,EUV 

SQb 

■  Ch(>|/ )  - 1]5(», V ),  b  =  e,p,  EUV 

5Qb\dz  J 


(A14) 

(A15) 
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b^e,p 


(A16) 


d  (dp, 


where 


b  =  e,p 


(A17) 


(A18) 


Explicitly 


-'O.b  n=l 


b  =  e,p 


^  (^)  1 

=  Pm^  y  0^  In  10  +(y  j"*  +  2y  log,o  Kt)-^  .  b  =  e,p 

L  ^0.6  J 


-'o.i  e 


9b 


b  =  e,p 


Now  we  produce  the  derivatives  with  respect  to  the  parameters 

dp,  ^  f  /) )  r  ~U.  ^  1  e~ti  « \  diif,  dh„ 


p(*). 

dp, 

P  n  * 

api"' ' 

or 

dp, 

ep" 

tr  ' 

-a^i^-'-llSK.O)  ,  b  =  e,p 


(A  19) 


(A20) 


(A21) 


^ ^=^’P  (A22) 


(A23) 


Similarly 
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an',*’  & 


ipr 

>  1 

1 _ 

b  =  e,p 


(A24) 


K: 


^  =  -[e'“”'Ch(x)-l]-^ 

oh,  '■  ■' 


(BVV) 


H 


Ch’(x)-3e"'‘'  Ch{x)  +  l];<‘''‘^ 


(A25) 

(A26) 


H- 


ap,  _  M, 


EUV 


dH_  H. 


(A27) 


II 

' 

loglO 

. 

1 - 1 

du.  Vi.  dHi,  Hh 

and  since  — -  =  — —  and  — 

BH,  H,  Bp':'  p':' 

-,  we  have 

’•  'KK’ 

f 

'  loglO 

V 

InlO 


■^0,6 


InlO 


Similarly 


[M,e"'"‘+(l-3?/,K"‘-l+w,] 


,(*) 


log 


^O.fr 


(A30) 


(A31) 


InlO,  b=e,p  (A32) 


dp,  _  (■£!),/,) 


b  =  e,p 


(A33) 
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b  =  e,p 


(A34) 


with  =  1  keV  and  =  8  keV. 


pr^(O): 


ap,.  _pr>(x)  . 

a  (dp^\  ^  1  di^r' 

dr^\Q){dzr r^r\<^)  dz 


(A35) 

(A36) 


P*: 


^Pfc  5pfc 


5p  A  i 


//,  a//,  ap,  h,  p,  ^ 


-l]5(i/j,0),  b  =  e,p 
-3e-’'‘+l)?(w„0),  b 


e,P 


(A37) 

(A3  8) 


While 


^P,  _  ^ 

attn  aa 


0  V 


dz 


=  0  and 


both  /;„£  and  are  explicitly 


dependent  on  a(r)  so 


a/;„,£  __  dh„,E  _  ^ 
az„  da. 


(A39) 


a//. 


0,  z  <  z„ 


'd^P 

^  1  ^P 

dz^ 

^=l'mE  _ 

(A40) 


dK,E  p, 

dz^  a/4 

dK,E  N„,E 

da,  7/x, 


dN„,E 


^  K.E-Z, 

2  Hi 


NE 


(A41) 

(A42) 

(A43) 
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